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1, 


SUMMARY 


The  main  goal  of  this  project  was  to  develop  an  efficient  clock  synchronization  scheme  to  ensure 
robust  operation  of  wireless  airborne  networks  in  conditions  of  arbitrary  network  delays  and 
absence  of  GPS  (Global  Positioning  Systems).  To  cope  with  the  Gaussian  or  non-Gaussian  nature 
of  the  random  network  delays,  a  novel  method,  referred  to  as  the  Gaussian  Mixture  Kalman 
Particle  Filter  (GMKPF),  is  proposed  herein  to  estimate  the  clock  offset  in  wireless  sensor 
networks.  GMKPF  represents  a  better  and  more  flexible  alternative  to  the  Gaussian  Maximum 
Likelihood  (GML),  and  Exponential  Maximum  Likelihood  (EML)  estimators  for  clock  offset 
estimation  in  non-Gaussian  or  non-exponential  random  delay  models.  The  computer  simulations 
illustrate  that  GMKPF  yields  much  more  accurate  results  relative  to  GML  and  EML  when  the 
network  delays  are  modeled  in  terms  of  a  single  non-Gaussian/non-exponential  distribution  or  as 
a  mixture  of  several  distributions.  As  deliverables,  the  set  of  Matlab  programs  used  to  implement, 
simulate  and  validate  the  performance  of  GMKPF  together  with  simulation  results  are  provided. 
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2. 


INTRODUCTION 


Future  wireless  airborne  networks  are  envisioned  to  represent  the  next  frontier  of  networking,  to 
be  pervasive  and  ubiquitous,  and  to  provide  a  wide  range  of  services  and  applications.  This  trend 
is  underlined  by  a  number  of  technological  advances  and  demands.  The  rapidly  growing  demands 
for  mobility  and  anywhere-anytime  data  access  represent  a  major  driving  force  behind  the  next 
generation  of  mobile  wireless  airborne  networks.  Recent  technological  developments  mark  also 
the  departure  of  telecommunications  systems  from  homogeneous  networks  to  heterogeneous 
networks,  from  non-intelligent  devices  to  smart  devices,  and  from  telephony-based  services  to 
multi-media  services.  In  addition,  recent  advances  in  hardware  and  inexpensive  wireless  radio 
systems  have  made  also  possible  the  design  of  low-cost,  low-power,  and  multi-functional  sensor 
devices.  When  deployed  in  a  large  number  across  a  geographical  area,  these  sensor  devices  create 
a  self-organized  cooperative  ad-hoc  network  that  is  perfectly  fit  for  distributed  sensing  and 
automated  information  gathering,  processing  and  communication.  The  upcoming  years  will  very 
likely  witness  a  growing  demand  for  more  intelligent  sensor  systems  that  will  be  networked  with 
wireless  local  area  networks  (WLANs),  Internet,  satellite  and  Unmanned  Aerial  Vehicle  (UAV) 
networks  to  create  a  global  wireless  airborne  network  with  increased  functionality  and 
performance. 

In  general,  for  distributed  computing  and  networking  systems,  maintaining  the  logical  clocks  of 
the  computers  in  such  a  way  that  they  are  never  too  far  apart  is  one  of  the  most  complex 
problems  of  computer  engineering.  Whether  it  is  the  disciplining  of  computer  clocks  with  the 
devices  synchronized  to  a  Global  Positioning  System  (GPS)  satellite  or  a  Network  Time  Protocol 
(NTP)  time  server  over  the  Internet,  it  is  possible  to  equip  some  primary  time  servers  for  the 
purpose  of  synchronizing  a  much  larger  number  of  secondary  servers  and  clients  connected 
through  a  common  infrastructure.  In  order  to  do  this,  a  distributed  network  clock  synchronization 
protocol  is  required  through  which  a  server  clock  can  be  read,  the  readings  to  other  clients  can  be 
transmitted  and  each  client  clock  can  be  adjusted  as  required.  In  such  a  distributed 
synchronization  approach,  the  participating  devices  exchange  timing  information  with  their 
chosen  reference  at  regular  intervals  and  adjust  their  logical  clocks  accordingly. 
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A  computer  clock  in  general  has  two  components,  namely  a  frequency  source  and  a  means  of 
accumulating  timing  events  (consisting  of  a  clock  interrupt  mechanism  and  a  counter 
implemented  in  software).  The  implementation  of  the  computer  clock  in  the  operating  system  and 
the  programming  interface  differ  between  operating  systems  and  hardware  platforms.  However, 
the  basic  source  of  timing  is  an  uncompensated  quartz  crystal  oscillator  and  the  clock  interrupts  it 
generates.  Theoretically,  two  clocks  would  remain  synchronized  if  their  offsets  are  set  equal  and 
their  frequency  sources  run  at  the  same  rate.  However,  practical  clocks  are  set  with  limited 
precision  and  the  frequency  sources  run  at  slightly  different  rates.  In  addition,  the  frequency  of  a 
crystal  oscillator  varies  due  to  initial  manufacturing  tolerance,  aging,  temperature,  pressure  and 
other  factors.  Because  of  these  inherent  instabilities,  distributed  clocks  must  regularly  be 
synchronized  to  keep  them  running  close  to  each  other. 

The  Network  Time  Protocol  [22],  [23]  represents  the  most  widely  used  clock  synchronization 
protocol  for  large-scale  networks  with  static  topology  such  as  the  Internet.  In  NTP,  the  nodes  are 
externally  synchronized  to  a  global  reference  time  that  is  represented  in  the  network  by  a  set  of 
master  nodes  or  time  servers  that  are  referred  to  as  layer- 1  servers.  The  entire  synchronization 
process  assumes  a  hierarchical  tree  organization  of  the  network  nodes.  Despite  its  wide-spread 
use  in  the  synchronization  of  Internet,  NTP  is  not  appropriate  for  synchronization  of  wireless  ad- 
hoc  sensor  networks  that  are  subject  to  severe  energy-constraints,  dynamic  topologies  caused  by 
mobility  and  node  failures,  and  absence  of  GPS  and  global  time  references  (due  to  either 
jamming,  interferences,  or  absence  of  direct  line  of  sight  communication  links).  In  addition,  the 
service  provided  by  NTP  assumes  continuous  synchronization  of  all  the  network  nodes  with 
maximum  accuracy  and  with  no  concern  about  energy  consumption.  However,  NTP  is  not 
equipped  with  a  mechanism  to  enable  the  local  synchronization  of  a  subset  of  nodes,  and  to  keep 
the  rest  of  the  nodes  switched  to  a  power-  saving  (sleeping)  state.  Since  listening  continuously  for 
the  synchronization  beacons  is  an  energy-consuming  operation,  NTP  cannot  directly  be  applied  to 
synchronization  of  energy-constrained  wireless  ad-hoc  networks  as  is  the  case  with  wireless 
airborne  networks  and  wireless  sensor  networks. 

These  considerations  illustrate  the  need  for  novel  distributed  and  scalable  synchronization 
protocols  for  wireless  ad-hoc  networks  that  in  general  must  satisfy  a  series  of  requirements: 
energy-efficiency,  robustness  with  respect  to  node  mobility  and  link/node  failures,  and  ability  to 
guarantee  the  long-term  network  synchronization  at  local  and  global  scales. 
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Clock  synchronization  is  important  for  many  applications  such  as  Internet  delay  measurements, 
cellular  networks,  data  security  algorithms.  Media  Access  Control  (MAC)  protocols  like  Time 
Division  Multiple  Access  (TDMA),  Internet  Protocol  (IP)  telephony,  ordering  of  updates  in 
database  systems,  etc.  Recently,  with  the  advent  of  Wireless  Sensor  Networks  (WSNs)  and 
Wireless  Airborne  Networks  (WANs),  developing  clock  synchronization  protocols  that  suit  their 
specific  requirements  is  becoming  an  important  research  problem.  A  large  number  of  their 
applications  require  the  clocks  of  the  nodes  to  run  synchronously  on  a  common  timescale.  This  is 
the  case  with  applications  such  as  data  fusion,  efficient  duty  cycling  operations,  acoustic 
beamforming,  localization,  security  and  object  tracking.  Unlike  conventional  networks,  energy 
efficiency  must  also  be  taken  into  account  for  addressing  the  clock  synchronization  problem  in 
WSNs  and  WANs. 
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3. 


METHODS,  ASSUMPTIONS  AND  PROCEDURES 


During  the  last  two  decades,  many  clock  synchronization  protocols  have  been  proposed  such  as 
[1],  [2],  [23],  etc.  NTP  [23]  is  a  protocol  for  synchronizing  the  clocks  of  computer  systems  over 
packet-switched,  variable-  latency  data  networks  and  it  represents  the  Internet  standard  for  time 
synchronization.  It  is  a  layered  client-  server  architecture  based  on  the  User  Data  Protocol  (UDP) 
message  passing  which  synchronizes  computer  clocks  in  a  hierarchical  way  using  the  offset  delay 
estimation  method.  NTP’s  sender-receiver  synchronization  architecture  is  widely  accepted  in 
designing  time  synchronization  algorithms  and  consists  of  the  same  two-way  timing  message 
exchange  mechanism  targeted  in  this  project. 

A  protocol  based  on  the  remote  clock  reading  method  was  put  forward  by  [2],  which  handles 
unbounded  message  delays  between  processes.  In  [1],  the  time  transmission  protocol  is  used  by  a 
node  to  communicate  the  time  on  its  clock  to  a  target  node,  which  subsequently  estimates  the 
time  in  the  source  node  by  using  message  timestamps  and  message  delay  statistics.  For  ad-hoc 
communication  networks,  the  time  synchronization  protocol  [8]  represented  one  of  the  pioneering 
contributions  in  this  area.  The  protocol  is  based  on  generating  timestamps  to  record  the  time  at 
which  an  event  of  interest  occurred.  The  timestamps  are  updated  by  each  node  using  its  local 
clock  and  the  time  transformation  method,  where  the  final  timestamp  is  expressed  in  terms  of  an 
interval  with  a  lower  bound  and  an  upper  bound.  In  the  realm  of  wireless  sensor  networks,  the 
clock  synchronization  protocols  of  particular  note  are  Reference  Broadcast  Synchronization  (RBS 
[5]),  Timing  Synch  Protocol  for  Sensor  Networks  (TPSN  [6])  and  Time  Diffusion  Protocol  (TDP 
[7]).  RBS  relies  on  simultaneous  reception  of  broadcast  pulses  by  several  nodes  transmitted  by  a 
common  neighboring  node  after  which  the  nodes  exchange  their  timestamps  and  estimate  the 
relative  time  offsets  and  skews.  On  the  other  hand,  TPSN  is  based  on  the  same  sender-receiver 
paradigm  as  in  NTP,  like  many  other  traditional  clock  synchronization  protocols.  The  basic 
difference  is  that  TPSN  has  been  molded  sufficiently  to  suit  the  requirements  of  wireless  sensor 
networks.  On  the  other  side,  TDP  establishes  a  network-wide  equilibrium  time  through  an 
iterative,  weighted  averaging  technique  based  on  a  diffusion  of  messages  involving  all  the  nodes 
in  the  synchronization  process. 
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Clock  synchronization  between  any  two  nodes  is  generally  accomplished  by  message  exchanges. 
Due  to  the  presence  of  non-deterministic  and  possible  unbounded  message  delays,  messages  can 
get  delayed  arbitrarily,  which  makes  the  clock  synchronization  very  difficult  [10].  The  most 
commonly  proposed  non-deterministic  network  delay  distributions  are  the  Gaussian,  exponential. 
Gamma,  and  Weibull  probability  density  functions  (pdfs)  [9],  [20],  [25].  In  general,  it  is  difficult 
if  not  impossible  to  assess  which  distribution  model  may  be  fit  to  capture  the  network  delay 
distributions  in  a  given  wireless  sensor  network  (WSN).  This  is  due  to  the  fact  that  various  factors 
might  impact  differently  the  distribution  of  network  delays  [17],  [18].  The  Gaussian  pdf  [12]  and 
the  exponential  pdf  [9]  were  also  recently  proposed  to  model  the  network  delays  in  WSNs. 
Herein,  the  maximum  likelihood  (ML)  estimators  for  clock  offset  estimation  in  the  presence  of 
Gaussian  and  exponential  network  delay  distributions  will  be  referred  to  as  the  Gaussian  ML 
(GML)  and  exponential  ML  (EML),  respectively.  Reference  [24]  shows  that  GML  and  EML  are 
quite  sensitive  to  the  network  delay  distributions.  Therefore,  one  important  problem  that  rises  up 
is  to  design  clock  offset  estimation  schemes  that  are  robust  to  the  distribution  of  unknown 
network  delays. 

To  overcome  these  challenges,  in  this  project  a  novel  clock  offset  estimation  method,  referred  to 
as  the  Gaussian  Mixture  Kalman  Particle  Filter  (GMKPF),  is  proposed  and  thoroughly  tested. 
Extensive  computer  simulations  illustrate  GMKPF’ s  merits  of  being  robust  and  yielding  very 
accurate  clock  offset  estimates  in  the  presence  of  arbitrary  network  delay  distributions.  The  clock 
offset  estimation  framework  adopted  in  this  project  is  identical  with  the  two-way  message 
exchanges  between  two  nodes,  encountered  in  NTP  [22]  and  TPSN  protocol  [24].  GMKPF 
combines  the  importance  sampling  (IS)  based  measurement  update  step  with  a  KF  (Kalman  Filter) 
based  Gaussian  sum  filter  for  the  time-update  and  proposal  density  generation.  Since  GMKPF 
employs  new  observations  and  exploits  the  Expectation-Maximization  (EM)  algorithm  to  obtain 
the  Gaussian  Mixture  Model  (GMM),  GMKPF  is  expected  to  exhibit  better  estimation 
performance  when  compared  to  GMF  and  EMF  in  general  non-Gaussian/non-exponential  delay 
models.  Thus  far,  in  the  synchronization  literature  for  WSNs,  it  appears  that  only  very  few 
preliminary  and  straightforward  applications  of  standard  Kalman  filtering  or  general  adaptive 
signal  processing  techniques  were  reported  (see  [14],  [16]  and  [26])  to  improve  the  mean  square 
error  (MSE)  performance  of  protocols  such  as  RBS  [12]  or  TPSN  [15]. 
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In  this  project,  upon  designing  the  GMKPF,  a  thorough  performance  analysis  of  GMKPF,  GML, 
and  EML  in  the  presence  of  the  two-way  message  exchange  mechanism  between  two  nodes  and 
symmetric/asymmetric  Gaussian,  exponential.  Gamma,  Weibull  network  delay  distributions  is 
first  carried  out.  The  performance  of  GMKPF,  GML,  and  EML  is  also  simulated  under  the 
mixing  of  two  different  distributions:  Gaussian  and  exponential,  Gaussian  and  Gamma,  Gaussian 
and  Weibull,  exponential  and  Gamma,  exponential  and  Weibull,  and  Gamma  and  Weibull  delay 
distributions,  respectively.  The  computer  simulation  results  corroborate  the  superior  performance 
of  the  proposed  method  relative  to  GML  and  EML,  and  its  robustness  to  general  network  delay 
distributions.  Therefore,  the  proposed  GMKPF  method  represents  a  high-performance  and  very 
reliable  clock  offset  estimation  scheme  fit  to  overcome  the  uncertainties  caused  by  the  network 
delay  distributions. 


7 


4, 


RESULTS  AND  DISCUSSION 


4.1.  Problem  Formulation  and  Objectives 

The  two-way  timing  message  exchange  mechanism  is  a  recently  proposed  clock  synchronization 
scheme  for  wireless  sensor  networks  [15],  [24].  Under  this  mechanism,  the  synchronization  of 
two  nodes  A  and  B  is  achieved  through  a  number  of  N  cycles.  Each  cycle  assumes  two  message 
transmissions:  one  from  node  A  to  node  B,  followed  by  a  reverse  transmission  from  node  B  to 
node  A.  At  the  beginning  of  the  k\h  cycle,  the  node  A  sends  its  time  reading  7]  ^  to  Node  B, 

which  records  the  arrival  time  of  the  message  as  Tj  ^  ,  according  to  its  own  time  scale.  Similarly  a 
time  message  exchange  is  performed  from  Node  B  to  Node  A.  At  time  T.^  ^  node  B  transmits 
back  to  node  A  the  time  information  and  Tj  ^ .  Denoting  by  7],  the  arrival  time  at  node  A  of 
the  message  sent  by  node  B,  node  A  would  then  have  access  to  the  time  information  T.  ^,j  = 
1,  .  .  .  ,  4  at  the  end  of  the  ^h  cycle,  which  provide  sufficient  information  for  estimating  the  clock 
phase  offset  6^  of  node  A  relative  to  node  B  clock.  The  message  exchanges  that  take  place 
between  two  generic  nodes  A  and  B  are  depicted  in  Fig.  1. 


9^  ;  clock  ofE^t  (time  difference) 


nodeB 


node  A 


Figure  1.  (a)  A  message  exchange  between  two  nodes  A  and  B  that  present  only  clock  phase  offset,  (b)  Multiple 
message  exchanges  between  nodes  A  and  B  that  present  clock  phase  offset  and  skew. 
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Utilizing  the  derivation  presented  in  [24],  the  differences  between  the  ^h  up  and  down-link  delay 
observations  corresponding  to  the  timing  message  exchange  are  given  by 

U  k'-=T2k-T^k  =  d  +  6^+  Lk  and  =T^k~d'2k  =d  -6^  +  Mk ,  respectively.  The  fixed  value  d 
denotes  the  fixed  (deterministic)  propagation  delay  component  (which  in  general  is  neglected 
«  0  in  small  range  networks  that  assume  RF  transmissions).  Parameters  Lk  and  M k  stand  for 

the  variable  portions  of  the  network  delays,  and  are  assumed  to  be  any  distributions  such  as 
Gaussian,  exponential.  Gamma,  Weibull  or  mixtures  of  two  different  distributions. 

Given  the  observation  samples  z^.  =  our  goal  is  to  find  the  minimum  variance  estimate 

of  the  unknown  clock  offset  9^ .  For  convenience,  the  notation  :=  9^  will  be  used  henceforth. 
Thus,  it  turns  out  that  we  are  looking  to  determine  the  estimator: 

Xk=E{Xk\7J},  (1) 

where  Z!  denotes  the  set  of  observed  samples  up  to  time  I,  =  {Zq,Zj,---,Z;}  .  Since  the  clock 

offset  value  is  assumed  constant,  the  clock  offset  can  be  modeled  as  obeying  a  Gauss-Markov 
dynamic  channel  model  of  the  form: 


^k 


FXk-i+Vk-i, 


(2) 


where  F  is  the  state  transition  matrix  for  clock  offset.  The  additive  noise  component  can  be 
modeled  as  Gaussian  with  zero  mean  and  covariance  £'{v^vj}  =  Q  .  The  vector  observation 
model  follows  from  the  observed  samples  and  it  assumes  the  expression: 


d  X.  Lk 

"11 

d  + 

■  1  ■ 

= 

d-Xk+Mk^ 

1 

-1 

(3) 


where  the  observation  noise  vector  accounts  for  the  random  delays.  One  can 

now  observe  that  eqs.  (2)  and  (3)  recast  our  initial  clock  offset  estimation  problem  into  a  Gauss- 
Markov  estimation  problem  with  unknown  states. 
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4.2.  A  Composite  Particle  Filtering  Approach 


Particle  filtering  is  a  sequential  Monte  Carlo  sampling  method  built  within  the  Bayesian  paradigm. 
From  a  Bayesian  perspective,  at  time  k,  the  posterior  distribution  p{Xj^  \  Zg,^)  is  the  main  entity 
of  interest.  Flowever,  due  to  the  non-Gaussianity  of  the  model  (3),  the  analytical  expression  of 
p{Xi^  I  Zg.^)  cannot  be  obtained  in  closed- form  expression,  excepting  for  some  special  cases  like 
Gaussian  or  exponential  pdfs.  Alternatively,  particle  filtering  can  be  applied  to  approximate 
piXf^  I  Zg,^)  by  stochastic  samples  generated  using  a  sequential  importance  sampling  strategy. 

Since  the  particle  filtering  with  the  prior  importance  function  employs  no  information  from 
observations  in  proposing  new  samples,  its  use  is  often  ineffective  and  leads  to  poor  filtering 
performance.  Flerein,  we  implement  a  slightly  changed  version  of  the  Gaussian  Mixture  Sigma 
Point  Particle  Filter  (GMSPPF)  proposed  in  [21],  and  which  will  be  referred  to  as  a  composite 
approach.  This  composite  approach  comes  out  from  the  utilization  of  another  filtering  technique 
producing  a  filtering  probability  density  function  used  as  importance  function  (IF)  for  the  particle 
filtering. 

The  GMSPPF  is  a  family  of  methodologies  that  use  hybrid  sequential  Monte  Carlo  simulation 
and  a  Gaussian  sum  filter  to  efficiently  estimate  posterior  distributions  of  unknown  states  in  a 
non-linear  dynamic  system.  However,  in  our  state  space  modeling,  because  of  the  linear  model, 
we  do  modify  this  method  further.  Following  [21],  we  will  next  describe  briefly  the  general 
framework  assumed  by  the  GMKPF  method,  obtained  by  replacing  the  SPKF  with  a  KF.  We  next 
outline  the  main  features  of  the  proposed  approach.  First,  we  remark  that  any  probability  density 
can  be  approximated  as  closely  as  desired  by  a  Gaussian  mixture  model  (GMM)  of  the 
following  form  [11], 


P(x)  ~  Pg{x)  =  (4) 

g=i 

where  G  stands  for  the  number  of  mixing  components,  denote  the  mixing  weights  and 
N(x;//,P)  is  a  normal  distribution  with  mean  p  and  covariance  P  .  Thus,  the  predicted  and 

updated  Gaussian  components,  i.e.,  the  means  and  covariances  of  the  involved  probability 
densities  (posterior,  importance,  and  so  on)  are  calculated  using  the  Kalman  filter  (KF)  instead  of 
the  Sigma  Point  Kalman  Filter  (SPKF)  [19],  [21].  Since  the  state  and  observation  equations  are 
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linear,  the  KF  was  employed  instead  of  the  SPKF.  Therefore,  the  resulting  approach  is  called  the 
Gaussian  mixture  Kalman  particle  fdter  (GMKPF).  In  order  to  avoid  the  particle  depletion 
problem  in  cases  where  the  observation  (measurement)  likelihood  is  very  peaked,  the  GMKPF 
represents  the  posterior  density  by  a  GMM  which  is  recovered  from  the  re-sampled  equally 
weighted  particle  set  using  the  Expectation-Maximization  (EM)  algorithm. 


In  general  for  the  particle  fdtering  approach,  the  posterior  density  |  Zj.^)  ,  where 

Xq.^  =  {Xo,---,x^}  and  Zj,^  =  ,  constitutes  the  complete  solution  to  the  sequential 

estimation  problem.  Our  objective  is  to  generate  samples  from  the  distribution  \  Zj.^).  For 

this  purpose,  we  have  collected  N  sets  of  samples  x^l  =  with  weights 

i  =  The  particles  approximate  pix^./,  IZj-^)-  Finally,  the  conditional 

mean  state  and  the  corresponding  error  covariance  can  be  calculated: 


N  N 


(5) 


1=1 


1=1 


At  the  end  of  each  recursion,  the  particles  are  resampled  to  ensure  they  occur  with  the  same 
probability  as  the  weights. 

The  GMKPF  combines  the  importance  sampling  (IS)  based  measurement  update  step  with  a  KF 
based  Gaussian  sum  fdter  for  the  time-update  and  proposal  density  generation.  In  the  time  update 
stage,  GMKPF  approximates  the  prior,  proposal  and  posterior  density  function  as  GMMs  using 
banks  of  parallel  KFs.  The  updated  mean  and  covariance  of  each  mixand  follow  from  the  KF 
updates.  In  the  measurement  update  stage,  the  GMKPF  uses  a  finite  GMM  representation  of  the 
posterior  filtering  density 


Pgix,  \z,)  =  Y,ccfN{xp,pf,Pt^^),  (6) 

^=1 

where  G  is  the  number  of  GMMs,  are  the  mixing  weights  and  N {xp,  is  a  normal 

distribution  determined  from  the  gth  KF  with  predicted  mean  =  x^  and  positive  definite 

covariance  .  This  is  recovered  from  the  weighted  posterior  particle  set  of  the  IS  based 

measurement  update  stage,  by  means  of  an  Expectation-Maximization  (EM)  [13]  step.  The  EM 
algorithm  can  be  used  to  obtain  Gaussian  Mixture  approximations  from  these  particles  and 
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weights.  Through  this  mechanism,  the  EM-based  posterior  GMM  further  mitigates  the  “sample 
depletion”  problem  through  its  inherent  “kernel  smoothing”  nature.  The  EM  algorithm  provides 
an  iterative  method  to  estimate  8  via 

6  =  arg  max  p{ju\6) ,  (7) 

0 

with  the  Gaussian  mixture  specified  by  the  parameter  set 
0  =  {af’ .  Specifically,  the  EM  algorithm  is  a  two-step 

iterative  algorithm  which  works  as  follows:  given  a  ,  it  finds  the  next  value  via 

•  E-step:  Q{9  \  9^’'^)  =  £'{log  p{ju  \  9)  \  9^^^} 

•  M-step:  =  arg  max  Q{9  \  9^^^) 

e 

The  reader  is  directed  to  reference  [13]  for  more  detailed  explanations  of  the  EM  algorithm  for 
GMM.  Finally,  the  conditional  mean  state  estimate  and  the  corresponding  error  covariance  can  be 
calculated  as  follows: 

Tt = ^Pk=Yj^k\Pk"^+ipf-wPk^-Xky^  (8) 

g=l  g=l 

Below  we  provide  a  fairly  pseudo-code  for  a  GMKPF  algorithm  that  is  fit  for  estimating  clock 
offsets  in  non-Gaussian  delay  models. 

Algorithm 

( 1 )  At  time  k- 1 ,  initialize  the  densities 

•  The  posterior  density  is  approximated  by 

Pgi^k-x  I  ) 

^=1 

•  The  process  noise  density  is  approximated  by 

Pg  ^^k-i ) = z  py\^^^k-x ;  p^i  ^  ) 

/=i 

•  The  observation  noise  density  is  approximated  by 

^ir  =  Z  ’  ^nf  .  ) 

7=1 
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(2)  Pre -prediction  step 


•  Calculate  the  pre-predictive  state  density  using  KF,  Pgix,^  \ 

•  Calculate  the  pre-posterior  state  density  using  KF,  Pg{x^.  |  z^) 

(3)  Prediction  step 

•  the  predictive  state  density  using  GMM,  p^  {x^.  \  z^  j ) 

•  Calculate  the  posterior  state  density  using  GMM,  Pgix^  \  z^.) 

(4)  Observation  Update  step 

•  Draw  N  samples  ^ ;  /  =  1,  •  •  • ,  M}  from  the  importance  density  function, 
9iXk\^k)  =  Pgi^kl^k) 

•  Calculate  their  corresponding  importance  weights: 


pi^klzDPgizn^k-i) 


PgiZkU^k) 


N 

•  Normalize  the  weights: 

i=\ 


•  Approximate  the  state  posterior  distribution  using  the  EM-algorithm, 
Pgi^kl^k) 

(5)  Infer  the  conditional  mean  and  covariance: 

•  At  =  z =  Z ^k^ ^Xk^  ~ At )^zT  - ^k y 

l=l  l=l 


Or  equivalently,  upon  fitting  the  posterior  GMM,  calculate  the  variables  in  eq, 

(8). 
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4.3.  Implementation  Aspects  of  GMKPF  Algorithm 


GMKPF  can  be  viewed  as  an  efficient  tool  to  perform  probabilistic  inference,  i.e.,  to  estimate  the 
hidden  variables  (states  or  parameters)  of  a  system  in  an  optimal  and  consistent  fashion  given 
noisy  or  incomplete  observations.  Fig.  2  depicts  a  general  framework  for  probabilistic  inference. 


Observed 


Unobserved 


Figure  2.  Probabilistic  Inference 


A  block  diagram  of  GMKPF  is  represented  in  Fig.  3. 


«-• 


Figure  3.  Constitutive  Blocks  of  Gaussian  Kalman  Particle  Filter 

The  Matlab  functions  and  routines  that  have  been  used  for  implementing  the  GMKPF  algorithm 
are  depicted  in  Figs.  4  and  5. 
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Main  function 

Asymmetnc  Exponential  main_gmm_asym_expo  m  and  g$sm_a$ym_expo  m 
Asymmetnc  Gaussian  matf)_gmm_asym_9auss  m  and  9ssin_asyfn_gauss  m 
Symmetric  Exponential  main_gmm_expo  m  and  gssm^expo  m 
Symmetric  Gaussian  mam^gmm^gauss  m  and  gssm^gauaa  m 
Gamma  mam^gmm^gamma  m  and  gssm^gamma  m 
Weibul  main^gmm^wei  m  and  gssm^wei  m 

Common  function 

GMMFIT  Fil  a  Gauatian  midura  modal  (GMM)  with  M  componania  to  dataaat  X  uaing  ao  aipactation  maamiization  algorithm  (EM) 
GMMPROBABNJTY  Calculataa  any  of  tha  ralatad  (Ihtough  Bayat  lula)  probabilitiaa  of  a  Gautaian  Mixiuia  Modal  (gmmOS)  and  a 
gwen  dataset  X  'piobjype' is  a  stnng  indicating  which  of  thefourpiobabildy 


% 

% 

p(xic)  p(C) 

likelihood  prior 

% 

P(Cp()  z -  postenorz  " 

— 

% 

P(X) 

evidence 

% 

GMMINmAUZE  :  Initialises  Gaussian  mixture  model  (GMM)  from  data 
KMEANS  Trains  a  k  means  cluster  model  Use  the  EM  algorithms 

GMMSAMPLE  DrawN  samples  from  the  Gaussian  mixture  modal  (GMM)  descnbad  by  the  GMM  data  structure 
SRCDKF  Square  Root  Central  Difference  Kalman  Filter  (Sigma-Pomt  Kalman  Filler  vanant)  (In  linear  model,  almost  same  as  Kalman 
Filter) 


Figure  4.  List  of  MATLAB  functions  used  to  impiement  the  Gaussian  Mixture  Kalman  Particle  Filter 


Time  update  Measurement  update 


0  ' 
p^(,x  I  z,)  =  Ta;^>N(x  ; 


Figure  5.  Main  MATLAB  Functions  used  to  implement  the  Gaussian  Mixture  Kalman  Particle  Filter 
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4.4.  General  Simulation  Results 


In  this  section,  computer  simuiation  resuits  wiii  be  offered  to  assess  the  performance  of  GMKPF, 
GML  [24],  and  EML  [24]-approaches  for  estimating  the  ciock  offset  in  wireiess  sensor  networks. 
We  consider  a  totai  of  10  delay  models:  asymmetric  Gaussian,  exponential.  Gamma,  Weibull,  and 
mixtures  of  Gaussian  and  exponential,  Gaussian  and  Gamma,  Gaussian  and  Weibull,  exponential 
and  Gamma,  exponential  and  Weibull,  and  Gamma  and  Weibull  distributions.  The  reason  for  this 
study  is  to  illustrate  that  the  proposed  method  is  robust,  exhibits  superior  performance  and  can  be 
applied  to  deal  with  any  delay  distribution.  The  stationary  process  is  assumed  to  achieve  a 
given  constant  variance  Q  =  \e  —  A  .  The  number  of  particles  and  GMM  are  100  and  3, 
respectively. 

Figs.  6-9  show  the  MSB  (Mean  Square  Error)  of  the  estimators  assuming  that  the  random  delay 
models  are  asymmetric  Gaussian,  exponential.  Gamma,  Weibull  pdfs,  respectively.  The 
subscripts  1  and  2  are  used  to  differentiate  the  parameters  of  delay  distributions  corresponding  to 
uplink  and  downlink,  respectively.  As  an  example,  the  parameters  Cj  and  (J^  in  Fig.  6  denote 

the  standard  deviations  of  uplink  and  downlink  asymmetric  Gaussian  network  delay  densities, 
respectively.  The  MSEs  are  plotted  against  the  number  of  observations,  ranging  from  5  to  25. 
Note  that  the  GMKPF  performs  much  better  (a  reduction  of  MSB  with  over  100%)  when 
compared  to  GME  or  EME.  It  is  interesting  to  note  that  the  MSB  of  GME  exhibits  better 
performance  than  EME  in  the  asymmetric  Gaussian  delay  model  case  and  poorer  performance  in 
the  presence  of  asymmetric  exponential.  Gamma,  and  Weibull  delay  models.  The  reason  for  this 
is  that  Gamma  and  Weibull  delay  models  are  closer  to  the  exponential  distribution  than  the 
Gaussian  distribution. 
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Asymmetric  Gajssian  delay  rixidel 


Figure  6.  MSEs  of  clock  offset  estimators  for  asymmetric  Gaussian  random  delays  (oi=l,  02=4) 


Figure  7.  MSEs  of  clock  offset  estimators  for  asymmetric  Exponential  random  delays  (Xi=l,  A2=5) 


17 


Mean  Square  Errors  Mean  Square  Errors 


Gamma  delay  model 


. 
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Figure  8.  MSEs  of  clock  offset  estimators  for  Gamma  random  delays  (ai=2,  Pi=l) 
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Weibull  delay  model 
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Figure  9.  MSEs  of  clock  offset  estimators  for  Weibull  random  delays  (ai=2,  Pi=2)  and  (a2=6,  ^2=2) 
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To  further  quantify  the  robustness  of  the  estimators,  we  studied  the  performanee  of  the  GMKPF, 
GML,  and  EML  under  various  network  delay  eonditions,  where  the  random  delay  models  are 
mixtures  of  two  distributions.  For  examples,  in  Fig.  10,  we  mix  equally  a  Gaussian  with  an 
exponential  delay  model,  each  having  a  weight  of  50%.  This  means  that  if  10  observations  are 
received,  5  observations  are  Gaussian  and  the  remaining  5  samples  assume  an  exponential 
distribution.  From  Figs.  10-15,  we  observe  that  GMKPF  clearly  outperforms  the  GML  and  EML. 
In  these  cases,  the  GML  presents  better  performance  than  EML  if  the  network  delay  process  is 
closer  to  a  Gaussian.  Otherwise,  the  EML  exhibits  better  performance  than  the  GML,  while 
GMKPF  outperforms  both  the  GML  and  EML. 
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Mixing  of  Gaussian  and  Exponential  delay  model 
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Figure  10.  MSEs  of  clock  offset  estimators  for  mixing  of  a  Gaussian  (oi=l,  02=1)  and  an  Exponential  (Ai=l,  A2=5) 
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Mixing  of  Gaussian  and  Gamma  delay  model 
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Figure  11.  MSEs  of  clock  offset  estimators  for  mixing  a  Gaussian  (Oi=l,  02=1]  and  a  Gamma  (ai=2,  |3i=2) 


Figure  12.  MSEs  of  clock  offset  estimators  for  mixing  a  Gaussian  (oi=l,  02=1)  and  a  Weibull  random  delay  (ai=2, 

Pi=2  and  02=6,  (32=2) 
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Figure  13.  MSEs  of  clock  offset  estimators  for  mixing  an  Exponential  (Ai=l,  ^2=5)  a  Gamma  {cli=2,  pi=5  and  a2=2 

P2=2) 
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Figure  14.  MSEs  of  clock  offset  estimators  for  mixing  an  Exponential  (Ai=l,  A2=5)  and  a  Weibull  (ai=2,  pi=2  and  02=6 

p2=2) 
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Figure  15.  MSEs  of  clock  offset  estimators  for  mixing  a  Gamma  (ai=2,  Pi=5  and  02=2,  ^2=2)  and  a  Weibull  (ai=2, 

Pi=2  and  02=6,  ^2=2) 
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4. 5.  In-Depth  Assessment  of  GMKPF  Performance 


•  TPSN  -  clock  offset  model  (1):  (Asymmetric  Gaussian) 

GMKPF  conditions:  Initial  value:  exponential  ML  value,  3-component  GMM,  500 
particles,  process  noise  variance:  le-6,  Monte  Carlo  Simulations:  200. 

The  MSEs  of  clock  offset  estimators  in  the  presence  of  asymmetric  Gaussian  delays  are 
presented  in  Fig.  16. 
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Figure  16.  MSEs  of  clock  offset  estimators  for  asymmetric  Gaussian  deiays 
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•  TPSN  -  clock  offset  model  (2):  (Asymmetric  Exponential) 

GMKPF  conditions:  Initial  value:  exponential  ML  value,  3-component  GMM,  500 
particles,  process  noise  variance:  le-6.  Monte  Carlo  Simulations:  200. 

The  MSEs  of  clock  offset  estimators  in  the  presence  of  asymmetric  Exponential  delays 
are  presented  in  Fig.  17. 
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Figure  17.  MSEs  of  clock  offset  estimators  for  asymmetric  Exponential  delays 
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•  TPSN  -  clock  offset  model  (3):  (Gamma) 

GMKPF  conditions:  Initial  value:  exponential  ML  value,  3-component  GMM,  500 
particles,  process  noise  variance:  le-6.  Monte  Carlo  Simulations:  200. 

The  MSEs  of  clock  offset  estimators  in  the  presence  of  asymmetric  Gamma  delays  are 
presented  in  Fig.  18. 
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Figure  18.  MSEs  of  clock  offset  estimators  for  asymmetric  Gamma  delays 
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•  TPSN  -  clock  offset  model  (4):  (Weibull) 

GMKPF  conditions:  Initial  value:  exponential  ML  value,  3-component  GMM,  500 
particles,  process  noise  variance:  le-6.  Monte  Carlo  Simulations:  200. 

The  MSEs  of  clock  offset  estimators  in  the  presence  of  asymmetric  Weibull  delays  are 
presented  in  Fig.  19. 
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Figure  19.  MSEs  of  clock  offset  estimators  for  asymmetric  Weibull  delays. 
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•  TPSN  -  clock  offset  model  (1):  (Gamma) 

GMKPF  conditions:  Initial  value:  exponential  ML  value,  3,  5,  7-component  GMM,  500 
particles,  process  noise  variance:  le-6.  Monte  Carlo  Simulations:  300. 

The  MSEs  of  clock  offset  estimators  in  the  presence  of  asymmetric  Gamma  delays  are 
presented  in  Fig.  20. 
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Figure  20.  MSEs  of  clock  offset  estimators  for  asymmetric  Gamma  delays. 
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•  TPSN  -  clock  offset  model  (2):  (Gamma) 

The  total  number  of  message  exchanges  in  a  network  with  100  nodes  and  MSE=0.001  are 
as  follows: 

Ntpsn^eml  =  2iV(Z  - 1)  =  2  X  920  X  (1 00  - 1)  =  1 82 1 60 
Ntpsn^gmkpf  =  2N{L  - 1)  =  2  X  30  X  (1 00  - 1)  =  5940 

The  number  of  timing  messages  vs.  MSEs  of  clock  offset  estimators  in  the  presence  of 
asymmetric  Gamma  delays  are  presented  in  Fig.  21.  Monte  Carlo  Simulations:  400. 
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Figure  21.  Number  of  timing  messages  vs.  MSEs  of  ciock  offset  estimators  for  asymmetric  Gamma  deiays 
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•  TPSN  -  clock  offset  model  (1):  (Weibull) 

GMKPF  conditions:  Initial  value:  exponential  ML  value,  3,  6,  10-component  GMM,  500 
particles,  process  noise  variance:  le-6.  Monte  Carlo  Simulations:  300. 

The  MSEs  of  clock  offset  estimators  in  the  presence  of  asymmetric  Weibull  delays  are 
presented  in  Fig.  22. 
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Figure  22.  MSEs  of  clock  offset  estimators  for  asymmetric  Weibuil  deiays 
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•  TPSN  -  clock  offset  model  (2):  (Weibull) 

The  total  number  of  message  exchanges  in  a  network  with  100  nodes  and  MSE=0.001  are 
as  follows: 

Ntpsn^eml  =  - 1)  =  2  X  220  X  (1 00  - 1)  =  399960 

Ntpsn^omkpf  =  ^N{L  - 1)  =  2  X  30  X  (1 00  - 1)  =  5940 

The  number  of  timing  messages  vs.  MSEs  of  clock  offset  estimators  in  the  presence  of 
asymmetric  Weibull  delays  are  presented  in  Fig.  23.  Monte  Carlo  Simulations:  400. 
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Figure  23.  Number  of  timing  messages  vs.  MSEs  of  clock  offset  estimators  for  asymmetric  Weibull  delays 
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5.  CONCLUSIONS 


Time  synchronization  is  a  significant  component  in  the  deployment  of  wireless  ad-hoc  networks, 
and  a  number  of  fundamental  operations,  like  data  fusion,  power  management  and  transmission 
scheduling,  require  accurate  time  synchronization.  Since  the  conventional  NTP  time 
synchronization  protocol  for  the  Internet  cannot  be  directly  applied  to  wireless  sensor  networks 
and  wireless  airborne  networks,  a  number  of  synchronization  protocols  have  been  developed  in 
this  project  to  meet  the  unique  requirements  of  these  applications. 

This  project  developed  a  very  general  and  powerful  inference  method  for  estimating  the  clock 
offset  in  wireless  ad-hoc  networks.  The  benefits  of  the  proposed  synchronization  method  are  in 
terms  of  improved  performance  and  applicability  to  arbitrary  random  delay  models  such  as 
asymmetric  Gaussian,  asymmetric  exponential.  Gamma,  Weibull,  as  well  as  to  mixtures  of  these 
delay  models.  One  negative  aspect  is  the  fact  that  analytical  closed  form  expressions  do  not 
necessarily  exist  and  in  general  it  is  hard  to  derive  lower  bounds  in  the  presence  of  (unknown) 
non-Gaussian  distributions.  The  project  proposed  a  robust  estimator  based  on  the  GMKPF  that  is 
capable  of  estimating  the  clock  offset  in  arbitrary  delay  models,  a  result  which  might  present 
applications  in  numerous  wireless  sensor  networks  applications  with  tight  synchronization 
requirements  [17],  [18].  Computer  simulations  also  show  that  the  proposed  method  yields 
superior  performance. 

The  proposed  statistical  inference  mechanism  is  quite  general  and  can  be  practically  applied  to 
any  modeling  problem  that  can  be  phrased  in  terms  of  a  state-space  representation.  The  additive 
noise  can  assume  arbitrary  distributions,  and  the  observation  equation  can  be  linear  or  nonlinear. 
In  the  presence  of  a  nonlinear  observation  equation,  a  slightly  more  general  modeling  framework 
in  terms  of  particle  filters  might  have  to  be  adopted  to  handle  the  nonlinearities.  However,  in  our 
present  study,  the  clock  observation  equations  are  linear  equations;  therefore,  a  bank  of  Kalman 
filters  was  sufficient  to  efficiently  track  the  unknown  clock  offset  parameters.  We  would  like  to 
emphasize  that  the  proposed  statistical  inference  engine  could  be  applied  to  more  general 
applications  such  the  problem  of  joint  synchronization  and  localization.  In  this  scenario,  one  can 
perform  two  tasks,  namely  clock  synchronization  and  localization  of  a  target  using  the  same  set  of 
signals/data  samples.  Therefore,  one  could  expect  high  performance  and  very  fast  algorithms  to 
be  developed  using  the  proposed  statistical  inference  mechanism.  Furthermore,  the  results  of  this 
project  could  be  extended  further  to  address  the  problem  of  joint  synchronization,  localization  and 
tracking  of  a  moving  target,  or  for  assessing  the  topography  of  a  wireless  network. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


AO  =  Always  On 

Adaptive  Sync  =  Adaptive  Synchronization 

CRLB  =  Cramer-Rao  Lower  Bound 

FTSP  =  Flooding  Time  Synchronization  Protocol 

GPS  =  Global  Positioning  System 

ID  =  Identification 

iid  =  independent  and  identically  distributed 

IP  =  Internet  Protocol 

MAC  =  Medium  Access  Control 

ML  =  Maximum  Likelihood 

MLE  =  Maximum  Likelihood  Estimator 

MSc  =  Master  of  Science 

MSE  =  Mean  Square  Error 

NTP  =  Network  Time  Protocol 

PBS  =  Pairwise  Broadcast  Synchronization 

PDF  =  Probability  Density  Function 

PhD  =  Doctor  of  Philosophy 

RBS  =  Reference  Broadcast  Synchronization 

RF  =  Radio  Frequency 

ROS  =  Receiver  Only  Synchronization 

RRS  =  Receiver  Receiver  Synchronization 

RV  =  Random  Variable 

SI  =  Sensor  Initiated 


SRS  =  Sender  Receiver  Synchronization 
TDMA  =  Time  Division  Multiple  Access 
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TDP  =  Time  Diffusion  Protocol 


TPSN  =  Time  Protocol  for  Synchronization  of  Sensor  Networks 

UAV  =  Unmanned  Aerial  Vehicles 

UDP  =  Universal  Data  Protocol 

WAN  =  Wireless  Airborne  Network 

WLAN  =  Wireless  Local  Area  Network 

wrt  =  with  respect  to 

WSN  =  Wireless  Sensor  Network 
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